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Abstract 

First principle calculation of the QCD spectral functions (SPFs) based on the lattice QCD 
simulations is reviewed. Special emphasis is placed on the Bayesian inference theory and 
the Maximum Entropy Method (MEM), which is a useful tool to extract SPFs from the 
imaginary-time correlation functions numerically obtained by the Monte Carlo method. 
Three important aspects of MEM are (i) it does not require a priori assumptions or 
parametrizations of SPFs, (ii) for given data, a unique solution is obtained if it exists, 
and (iii) the statistical significance of the solution can be quantitatively analyzed. 

The ability of MEM is explicitly demonstrated by using mock data as well as lattice 
QCD data. When applied to lattice data, MEM correctly reproduces the low-energy 
resonances and shows the existence of high-energy continuum in hadronic correlation 
functions. This opens up various possibilities for studying hadronic properties in QCD 
beyond the conventional way of analyzing the lattice data. Future problems to be studied 
by MEM in lattice QCD are also summarized. 
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1 Introduction 



The hadronic spectral functions (SPFs) in quantum chromodynamics (QCD) play an 
essential role in understanding properties of hadrons as well as in probing the QCD vacuum 
structure. For example, the total cross section of the e"*" + e~ annihilation into hadrons 
can be expressed by the spectral function corresponding to the correlation function of the 
QCD electromagnetic current. Experimental data on the cross section actually support 
the following picture: asymptotically free quarks are relevant at high energies, while the 
quarks are confined inside hadronic resonances (such as p,uj,(j)) at low energies. Another 
example of the application of SPF is the QCD spectral sum rules, where the moments of 
SPF are related to the vacuum condensates of quarks and gluons through the operator 
product expansion 0. 

The spectral functions at finite temperature (T) and/or baryon chemical potential (fi) 
in QCD, which were originally studied in ref. 0], are also recognized as a key concept 
to understand the medium modification of hadrons (see, e.g., @, H, |, The physics 
motivation here is quite similar to the problems in condensed matter physics or in nuclear 
physics: how elementary modes of excitations change their characters as T and/or /i 
is raised. The enhanced low-mass e~^e~ pairs below the p-resonance and the suppressed 
high-mass pairs at J/ip and ?/^'-resonance observed in relativistic heavy ion collisions 

at CERN SPS ^ are typical examples which may indicate spectral changes of the qq 
system due to the effect of the surrounding environment (see recent reviews, [IC, 11]). 

Since the numerical simulation of QCD on a lattice is a first-principle method having 
remarkable successes in the study of "static" hadronic properties (masses, decay con- 



stants, • • •, etc.) |1^, 0, it is desirable to extend its power also to extracting hadronic 
SPFs. However, Monte Carlo simulations on the lattice have difficulties in accessing the 
"dynamical" quantities such as spectral functions and the real-time correlation functions. 
This is because measurements on the lattice can only be carried out at a finite set of 
discrete points in imaginary time. The analytic continuation from imaginary time to real 
time using limited and noisy lattice data is not well-defined and is classified as an ill-posed 
problem. This is the reason why studies to extract SPFs so far had to rely on specific 
ansatze for the spectral shape |14, 15 1. 

In this article, we review a new way of extracting the spectral functions from the 
lattice QCD data and explore its potentiality in detail. The method we shall discuss is 
the maximum entropy method (MEM), which was recently applied to the lattice QCD 
data by the present authors [|16] . 

In the maximum entropy method. Shannon's information entropy [|r^ and its gener- 
alization (which we call Shannon- Jaynes entropy throughout this article) play a key role. 
After the pioneering works on the applications of the information entropy to statistical 
mechanics by Jaynes ||18| (see also |jl9[) and to optical image reconstruction by Frieden 



MEM has been widely used in many branches of science [^. Applications include 
data analysis of quantum Monte Carlo simulations in condensed matter physics [^, ^ 
and in nuclear physics [^, image reconstruction in crystallography and astrophysics [^ , 
and so forth. In the context of QCD, MEM is a method to make statistical inference on 
the most probable SPF for given Monte Carlo data on the basis of the Bayesian probability 
theory [|25l . 



2 



In the maximum entropy method, a priori assumptions or parametrizations of the 
spectral functions need not to be made. Nevertheless, for given lattice data, a unique 
solution is obtained if it exists. Furthermore, one can carry out error analysis of the ob- 
tained SPF and evaluate the statistical significance of its structure. Our basic conclusion 
is that MEM works well for the current lattice QCD data and can reproduce low-lying 
resonances as well as high-energy continuum structure. This opens up various possibilities 
for studying hadronic properties in QCD beyond the conventional way of analyzing the 
data.Q 

This article is organized as follows. 
In Section 2, we shall give a general introduction to the structure of the SPFs in QCD. 
In Section 3, the MEM procedure will be discussed in detail. This section is written in a 
self-contained way so that the readers can easily start MEM analysis without referring to 
enormous articles that have been published in other areas. The uniqueness of the solution 
of MEM is proved explicitly in this Section. A subtle point related to the covariance 
matrix of the lattice data is also discussed. 

In Section 4, we present the results of the MEM analysis using the mock data. It is also 
discussed how the resultant SPF and its error are affected by the quality of the mock 
data. 

In Section 5, SPFs at T = in the mesonic channels are extracted by MEM using quenched 

Monte Carlo data obtained on 20'^ x 24 lattice with (3iat = 6.0. Detailed analyses of the 

spectral structure in the vector and pseudo-scalar channels are given. Part of the results 

presented in Section 4 and Section 5 have been previously reported in |T6| . 

In Section 6, future problems of MEM both from technical and physical point of view are 

summarized. 

In Appendix A (B), the statistical (axiomatic) derivation of the Shannon- Jaynes entropy 
is shown. In Appendix C, the singular value decomposition of matrices, which is crucial 
for numerical MEM procedure, is proved. 



We should mention here that there were at least two attempts aiming at a goal similar to this article 
but with different methods |Q. Realistic analyses of SPFs for QCD, however, have not yet been carried 
out in those methods. 
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2 Sum rules for the spectral function 



2.1 Real and imaginary time correlations 

For simplicity, we shall focus on SPFs at vanishing baryon chemical potential. We start 
by defining the following two-point correlation functions at finite T: 



D^^,{t, x) 



Drin'ir, x) 



Tr (r [J^(t,x) 4(0'0 
- i> , 



(27r)4 

Tr(T, [j,(r,x) ^(0,0 
, X — ^ f d k 



-H/T 



,-H/T 



/z 



/z 



(2.1) 



(2.2) 



(27r)^ 



Here and J^, are composite operators in the Heisenberg representation in real (imag- 
inary) time for (D). r] and rj' denote possible Lorentz and internal indices of the 
operators. Z is the partition function of the system defined as Z = Tie^^^^. Un is the 
Matsubara frequency defined by Un = 2n7iT when and J^, are the Grassmann even ( = 
bosonic) operators and Un = {2n + l)7rT when they are Grassmann odd ( = fermionic) op- 
erators. The retarded product R and the imaginary-time ordered product T,- are defined, 
respectively, as 



R[A{t)B{0)] = 9{t){A{t)B{0)TB{0)A{t)), 
Tr[A{T)B{0)] = e{T)A{T)B{0)±e{-T)B{0)A{T). 



(2.3) 
(2.4) 

The upper sign is for bosonic operators and the lower sign is for fermionic operators, which 
is used throughout this article. Eq.( |2.1|) and eq. (|2.2|) are called the retarded correlation 
and the Matsubara correlation, respectively. 

One can of course define two-point correlations for more general composite operators, 
but Jr^ and J^, are sufficient for the later discussions in this article. 



2.2 Definition of the spectral function 

The spectral function (SPF) is defined as the imaginary part of the Fourier transform of 
the retarded correlation (see, e.g., p7[]). 



Ar,r^>{ko,k) = - Im D^,{ko,k) 



TT 



(2.5) 



T-»0 



{2nr E "-y-{n\UO)\m){m\jl,mn)il T e-^^--'^)5\k^ - P^J (2.6) 
(2vr)^E(0|^.(0)|m)(m|4(0)|0)5^(A;^-P^), (2.7) 



where |n) is a hadronic state with 4-momentum and the normalization (m|n) = 8nm- 
P^„ is defined by P!^ — P^. D^{ko, k) is a Fourier transform of D^(t, x) with respect to 
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both time and space coordinates. In eq.( |2.7| ), the energy of the vacuum is renormahzed 
to zero. 

When rj = ■}]', the SPF has special properties. First of all, it is positive semi-definite 
for positive frequency: 

A^^{ko>0,k)>0. (2.8) 

Also, it has a symmetry under the change of variables: 

^w(-^o, -k) = TArinih, k) = ^^^^(fco, -k), (2.9) 

where we have assumed parity invariance of the system. 

As already mentioned in Section 1, A^^/ is in some cases directly related to the exper- 
imental observables. For example, consider = j^"^ and J^^, = with 

2_ 1 - 1_ 

jT = 3"^M« - 3^7mC^ - g^T/^s ■ ■ ■ , (2.10) 

which is the electromagnetic current in QCD. Then, the standard i?-ratio in the + e~ 
annihilation is related to the SPF at T = (eq.(p.7[)) as follows (see, e.g., 0): 

47r2 3 

R{s) = -—Y.A^,{s), (2.11) 

with s = /cq — k^. The dilepton production rate from hot matter at finite T, which is an 
observable in relativistic heavy-ion collisions, is written through eq.( p.6|) as follows (see, 
e.g., fig): 



where a is the electromagnetic fine structure constant. The mass of leptons is neglected 
in this formula for simplicity. Note that this equation is valid irrespective of the state of 
the system, i.e., whether it is in the hadronic phase or in the quark-gluon plasma. 

2.3 Dispersion relations 



Using the definition of the two point correlations eqs.( |2.1|j2!^ ) together with the gen- 



eral form of the spectral function eq.( ^.5D , one can derive the dispersion relation in the 
momentum space p] 



duj "^''y > - (subtraction), (2.13) 
-oo oj — kQ — le 

/+00 ^ /(a;,k) 

duj — — — (subtraction). (2-14) 
-oo UO — iuJn 



The integrals on the right hand side (r.h.s.) of eqs.( p.l3" , p.l4 ) do not always converge and 



need appropriate subtractions. This is because, in QCD, 

A^^/(co' oo, k) oc cj", (2.15) 

where n = dim[J,,] + dim[J,,/] — 4, with dim[J] being the canonical dimension of the 
operator J. For example, n = 2 for the mesonic correlation with J qq and n = 5 for 
baryonic correlation with J ~ qqq. 
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2.4 Sum rules 



From the dispersion relations eqs. (|2.13 , 2.14|) , one can derive several sum rules for the 
spectral functions. 



For example, for bosonic operators with rj' = rj, eq.( |2.13| ) is rewritten as 
Dj^^(A;0,k) = Hduj' - (subtraction). 



(2.16) 



In the deep-Euclidean region —>■ — oo with finite k, we make the operator product 
expansion of the left hand side (l.h.s.) of (|2.16| ) and subsequently apply the following 
Borel transformation on both sides of eq.( p.l6| ): 



B 



M 



lim 

-fcn, n — ► oo 



(-*o)' 



(2.17) 



Then one arrives at the Borel sum rule in the medium 



^2^raJ-2 C_(a, (M^) ) -^^^^^ ^^^^^ 

ra,m=0 



M2« 




77^ • (2-18) 



Here the r.h.s. of ( p.l8| ) is an asymptotic expansion in terms of t^q^jj^jM'^ ^ T'^ /M'^ and 
k^/M^ with being the Borel mass. Therefore, should be large enough for the 
expansion to make sense. Cnm are the dimensionless coefficients calculated perturbatively 
with the QCD running coupling constant ^^(M^) defined at the scale M^. {{02n{M'^))) is 
the thermal average of all possible local operators with canonical dimension 2n renormal- 
ized at the scale M^. Note that T enters only through the matrix elements {{02n{M'^))) . 

contains not only the Lorentz scalar operators but also Lorentz tensor operators be- 
cause of the existence of a preferred frame at finite T. In an appropriate Borel window 
{Mmin < M < Mmax)i ( p.l8|) givcs coustraiuts on the integrated spectral function. The 
in- medium sum rule in the form, eq.( p.l8| ), was first derived in It is a natural gen- 
eralization of the original QCD sum rule in the vacuum [0, |^. Recent applications of the 
in- medium QCD sum rules can be found in p5[ ]. 

Another form of sum rule, which is closely related to the main topic of this article, is 
derived from eq. ( |2.14| ) . Let us define the mixed representation of the Matsubara correlator. 



k) 



e ^-"^Z}^,,(^a;„,k) 



(2.19) 



du Ajjj^i(uj, k) 
Then, by using the identity. 



ijj — lUJr, 



-XT 



^ x-iujn 1 =F e 
one arrives at the sum rule 



Dr^r)' {r, k) 



+ 00 



1 =F e-^'^ 



(0 < r < (3), 



A^ri'i^^, k) duj 



(0 < r < p). 



(2.20) 



(2.21) 
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Eq.( p.2lD is always convergent and does not require subtraction as long as r 7^ 0. This 
is because A{uj,'k) in QCD has at most power-like behavior at large u (see eq.( p.l5| )). 
Owing to this, we always exclude the point r = in the following analysis. 

When r] = 1]' and J,, is a quark bilinear operator such as qq, q'j^q, ■ ■ ■ etc., one can 
further reduce the sum rule into a form similar to the Laplace transform: 



-{f3-T)iJ 

A{uj, k) du 



1 - e-/3'^ 

r+oo 

= / K{t,uj) A{uj,k) du {0<T<f3), (2.22) 
Jo 

where we have suppressed the index t] for simplicity. K is the kernel of the integral 
transform. Eq.( |2.22D is the basic formula, which we shall utilize later. From now on, we 
focus on SPFs "at rest" (k = 0) for simplicity and omit the label k. 

To get a rough idea about the structure of D{t), let us consider a parametrized form 
of SPF, ^4^(0;), for the correlation of isospin 1 vector currents. It consists of a pole (such 
as the p-meson) + continuum: 

^2 c-/, ,2 2 ^ 1 - '^f 



p^u) = 2F^6iu'-mi) + — [l + -^)eiuj-^o). (2.23) 

Here and are the mass of the vector meson and the continuum threshold, respec- 
tively, is the residue at the vector meson pole. This simple parametrization of SPF has 
been commonly used in the QCD sum rules [jl], (For a more realistic parametrization, 
see eq. (|4.5| ) in Section 4.2.) D{t) at T = corresponding to (|2.23|) reads 

D(r) = F^^m^e-vr + (l + f ) (:| + + ^) e'-^ (2.24) 

At long distances (large r), the exponential decrease of (|2.24|) is dominated by the pole 
contribution. On the other hand, at short distance (small r), the power behavior 
from the continuum dominates eq.( |2.24| ). Although ( |2.24| ) captures some essential parts. 



D{t) in the real world or on the lattice has richer structure. Studying this without any 
ansatze like ( |2.23D is a whole aim of this article. 

2.5 An ill-posed problem 

Monte Carlo simulation provides D{r) in (|2.22| ) for a discrete set of points, t = Ti, with 



1 < n/a < Nr, (2.25) 

where Nr is the number of the temporal lattice sites and a is the lattice spacing. In the 
actual analysis, we use data points in a limited domain r G [Tmin,Tmax]- Since we can 
generate only finite number of gauge configurations numerically, the lattice data D{Ti) 
has a statistical error. From such finite number of data with noise, we need to reconstruct 
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the continuous function A{uj) on the r.h.s. of eg. ( p.22| ), or equivalently to perform the 
inverse Laplace transform. 

This is a typical ill-posed problem, where the number of data points is much smaller 
than the number of degrees of freedom to be reconstructed. The standard likelihood 
analysis (x^-fitting) is obviously inapplicable here, since many degenerate solutions appear 
in the process of minimizing x^- This is the reason why the previous analyses of the 
spectral functions have been done only under strong assumptions on the spectral shape 
[0, |T5|. Drawbacks of the previous approaches are twofold: (i) a priori assumptions on 



SPF prevent us from studying the fine structures of SPF, and (ii) the result does not have 
good stability against the change of the number of parameters used to characterize SPF. 
Both disadvantages become even more serious at finite T, where we have almost no prior 
knowledge on the spectral shape. 

The maximum entropy method, which we shall discuss in the next section, is a method 
to circumvent these difficulties by making a statistical inference of the most probable SPF 
(or sometimes called the image in the following) as well as its reliability on the basis of a 
limited number of noisy data. 



8 



3 Maximum entropy method (MEM) 



In this section, we shall discuss the MEM procedure in some detail to show its basic 
principle as well as to show crucial points in its application to lattice QCD data. 

The theoretical basis of the maximum entropy method is Bayes' theorem in probability 



theory ||25| 



PIA-IVI ^ . (3.1) 

where P[X|y] is the conditional probability of X given Y. The theorem is easily proved 
by using the product formula P[XY] = P[X\Y]P[Y] = P[Y\X]P[X]. Let D stand for 
Monte Carlo data with errors for a specific channel on the lattice and H summarize all 
the definitions and prior knowledge such as A{uj > 0) > 0. From Bayes' theorem, the 
conditional probability of having A{uj) given the data reads 

PiAlnm ^ ^Mp. (3,2) 

Here P[D\AH] and P[A|/7] are called the likelihood function and the prior probability, 
respectively. P[D\H] is simply a normalization constant independent of A{uj). (Note here 
that it may be more appropriate to call P "plausibility" instead of "probability" , since P 



does not necessarily have the frequency interpretation [IS 



Now, the most probable image is A(uj) that satisfies the condition, 

6P[A\DH] , ^ 

^Jj-J = 0. (3.3) 

Furthermore, the reliability of the image satisfying ( |3.3| ) can be estimated by the second 
variation, or schematically, b'^P\A\Dli\j bAbA. 

When more data become available, P\A\Dli\ can be updated. This is seen from the 
following chain rule, which is a consequence of Bayes' theorem and the product formula 
for conditional probabilities: 

P\A\DrD<,H\ = P[D2\DiAH]P[Di\AH]P[A\H]/P[DiD2\H]. (3.4) 



For further discussions on the general use of Bayesian analysis, see ref. . 

To go further, we need to specify the explicit forms of the likelihood function and the 
prior probability. This will be discussed below. 

3.1 Likelihood function 

For large number of Monte Carlo measurements of a correlation function, the data is 
expected to obey the Gaussian distribution according to the central limit theorem: 

P[D\AH] = ^e-^ , (3.5) 
L = \Y.mn)-DA{n))Cr^{D{T,)-DA{T,)), (3.6) 
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where i and j run over the actual data points which we utihze in the analysis, Tmin/a < 
hj < Tmax/o.. FoT later purposes, we define the number of data points to be used in 
MEM, 



D{Ti) is the lattice data averaged over gauge configurations, 

^ AT 



Din 



N, 



conf m=l 



(3.7) 



(3.8) 



where Nconf is the total number of gauge configurations and D^iji) is the data for the 
m-th gauge configuration. D^(rj) in ( |3.(j| ) is the correlation function defined by the r.h.s. 
of eq.(p:22D. 

C in (|3.6|) is an X iV covariance matrix defined by 



N, 



conf 



con fi^-^ conf 1) m=l 



(3.9) 



Lattice data have generally strong correlations among different r's, and it is essential to 
take into account the off-diagonal components of C. 

The integration of P[D\AH] over D with the measure [dD] defined below is normalized 
to be unity. Zl denotes the corresponding normalization constant: 



[dD] = n dD{n), Zl = i2TTf/'^y/dM. 



(3.10) 



In the case where P[y4|iJ] = constant in eg. ( p.2| ), maximizing P[A\DH] is equivalent 
to maximizing eq. ( p.5|) with respect to A, which is nothing but the standard x^-fitting. 
On the lattice, as we shall see later, the number of lattice data is (9(10), which is much 
smaller than the number of points of the spectral function to be reproduced ((9(10^)). 
Therefore, the x^-fitting does not work. This difficulty is overcome in the maximum 
entropy method, where the existence of a non-constant P[yl|if] plays an essential role. 



3.2 Prior probability 

In MEM, the prior probability is written with auxiliary parameters a and m as 
1 



aS 



P[A\Ham] = — e 
Zs 

where S is the Shannon- Jaynes entropy. 



S 



A{uj) - m(cj) - A{uj) loj 



Aiu) 



miuj] 



E 

1=1 



Ai 

Ai-mi- Ai log ( — 

m 



duj 



(3.11) 

(3.12) 
(3.13) 
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a is a real and positive parameter, while m(a;) is a real and positive function called the 
default model or the prior estimate. Although a and m are a part of the hypothesis 
H in ( |3.2| ), we write them explicitly on the l.h.s. of ( |3.11D to separate them from the 
other hypotheses. (To be consistent with this notation, we replace P[D\AH] in ( p.5| ) by 
P[D\AHam] in the following, although P[D\AHam\ does not depend on a or m.) a 
will be integrated out later and is eliminated in the final results, m remains in the final 
results, but one can study the sensitivity of the results against the change of m. 

In the numerical analysis, the frequency is discretized into A^^^ pixels of an equal size 
Au as shown in (|3.13| ); namely, Ai = A{u!i)Au! and mi = m{uJi)Auj with ui = I ■ Au 
il<l<N^). 

The integration of P[A\Ham] over A with the measure [dA] defined below is normal- 
ized to be unity. Zs is the corresponding normalization factor: 

dA, /27r\^"/2 

'"^i^n^. ^.-(1) . (3.14) 

In Appendix A, we give a derivation of the Shannon- Jaynes entropy S, the measure 
[dA] and the normalization Zs, on the basis of the so-called "monkey argument" [31, 
|3^ , |33| . The prior probability with eq. ( |3.12| ) is shown to be the most unbiased one for 
positive images. The structure of S in (|3.12|) can be alternatively derived on an axiomatic 



basis (see and references therein): For completeness, we present, in Appendix B, the 
axioms somewhat simplified from those given in p4| . 



3.3 Outline of the MEM procedure 



The procedure of the maximum entropy method may be classified into three classes: 
historic, classical and Bryan's method [^. They are different in the treatment of a 
as well as the way to search the maximum of P[A\DHam\. In this article, we shall follow 
Bryan's method, which is the state-of-art MEM with the most efficient algorithm and the 
least conceptual difficulty. The method consists of the steps given below. 

Step 1: Searching for the most probable image for a given a. 
Combining 



and ( 3.11 ), one obtains 



P[A\DHam] oc 



Q{A) = aS-L. 



(3.15) 



Therefore, the most probable image for a given a (and m), which we call A^,, satisfies 



SQ 



SA{uj) 



0. 



(3.16) 



A — Aq 



At this stage, a plays a role of a parameter which controls the relative weight of the 
entropy S (which tends to fit A to the default model m) and the likelihood function L 
(which tends to fit A to the lattice data). 

As will be proved in Section 3.5, the solution of eq.( 3.16| ) is unique if it exists. This 
makes the MEM analysis robust and essentially different from the x^-fitting. The latter 
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can have many degenerate solutions in ill-posed problems. Further details on the algorithm 
for solving (|3.16|) are given in Section 3.4. 

Step 2: Averaging over a. 

The final output image Aout is defined by a weighted average over A and a: 

Aout{^) = J [dA] J da A{uj) P[A\DHam\P[a\DHm\ 

da Aa{uj) P[a\DHm] . (3.17) 



Here, we have assumed that P[A\DHam\ is sharply peaked around Aa{uj), which should 
be satisfied for good data, i.e., data with small errors. Under this assumption P[a\DHm] 
can be evaluated using Bayes' theorem as 

P [a I DHm] = J [dA] P [D \ AH am] P [A \ Ham] P [a \ Hm] /P [D \ Hm] 

oc P[a\Hm] [[dA] -^e^^^^ (3.18) 
J ZsZl 



oc P[a\Hm] exp 



1 a 
Y^log^- + aS{A^)-L{A^) 



(3.19) 



Here, A^'s are the eigenvalues of the real symmetric matrix in the frequency space. 



dA,dA,, 



(3.20) 



A=Aa 



The standard choice of the prior probability for a is either Laplace's rule {P[a\Hm] = 
const.) or Jeffreys' rule {P[a\Hm] = l/a) However, the integral on the r.h.s. of 

eq. ( |3.17| ) is insensitive to the choice, as long as the probability is concentrated around its 
maximum at a = a. We have checked that this is indeed the case for our lattice QCD 
data. Therefore, we use the Laplace rule for simplicity throughout this article. 

As for the averaging over a, we first determine a region of a, [amin,o:max], by the 
criterion P[a\DHm] > 10~^ x P[a\DHm]. Then, after renormalizing P[a\DHm] so that 
■fa^in P[a\DHm] = 1 is satisfied, we carry out the integration eq.( |3.17|) over the above 
interval. 

In the analysis in Section 4 and Section 5, we make a further approximation on 
P[a\DHwi\: After obtaining a by maximizing eq. ( 3.19|) , Aa and \k{oi) in ( ^.191) are 



assumed to have weak a-dependence and are replaced by A^ and \k{ci) respectively. The 
approximate form P"PP[a|Difm] is used for carrying out the averaging in ( p. 171) . We have 
checked that this approximation does not lead to an error larger than the line-width of 
the figures for SPF, although P°-'P'P[a\DHm\ itself for large a differs from P[a\DHm] by 
as much as 10%. 

Step 3: Error analysis. 

The advantage of MEM is that it enables one to study the statistical significance of 
the reconstructed image Aioj) quantitatively. The error should be calculated for A{ijj) 
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averaged over some interval in a;, since there are correlations among A{uj) at neighboring 
cj's. This is explicitly seen from the Hesse matrix of Q defined by 



6A{iu)6A{uj') 



(3.21) 

A.=A.cv 



which is not diagonal in general. 

In the following MEM analysis, the error estimate is carried out for (weighted) averages 
of the image Aout{^) in finite regions instead of estimating errors at each pixel For 
this purpose, we first define the average of A^u) given a, 

_ J[dA] J.du; A{u;)P[A\DHam]Wicu) J.dco A^{cu)W{uj) 
^^"^^''^ = J.du; Wiu) ^ JjdcoWiu;) ' ^^"^^^ 

where / = [co'i,ti;2] is a given region in the w-space and W{uj) is a weight function. To 
get the last expression in eg. ( |3.22| ), it is assumed as before that PlAlDHam] is highly 
peaked around Aa{uj). 

For simplicity we use the flat weight, W{lj) = 1, and omit the suffix W from now on. 
The variance of {Aa)i is similarly estimated as 

{{SA^f)i = f [dA] f duduj' 6A{uj)6A{uj')P[A\DHam] / f dudu' (3.23) 

d^dcu' (ttT^^TTTI^I I [ diuduj', (3.24) 

ixi \6A{uj)6A{uj') J ^^^^ Jixi 

where 6A{u) = A{u) — Aa{uj). The Gaussian approximation of P[A\DHam] around A^ 
is taken in the last expression. 

As we have done for Aout{^), the error for Aout in the region / is also given by the 
following average; 

{iSAout)^)i = j da {{5A^f)i P[a\DHm] . (3.25) 

By using the same procedure with a slight modification, it is possible to estimate 
the errors for quantities derived from the reconstructed image. One example is the pole 
residue, j^^i^ Aoutil^) ^^^^ with ''pole" implying the integral over the lowest isolated pole. 



Step 4: Sensitivity test under the variation of m. 

One can study the sensitivity of the final results under the variation of m. In QCD, 
m{uj large) can be estimated with perturbation theory (see eq. (|5.16|) ). Therefore, the 
value of m is known at least at large uj. 

When the final result is not stable enough against the variation of m, one may select 
the optimal m that gives the image with the smallest error for the spectral functions. 
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3.4 Maximum search using SVD 

The non-trivial part of the MEM analysis is to find the global maximum of Q in the 
functional space of A{uj), which has typically O(IO^) degrees of freedom in our case. 
Fortunately, the singular value decomposition (SVD) of the kernel K reduces the search 
direction into a subspace with a dimension not more than the number of data points 
~ 0(10). This has been shown by Bryan In the following, we shall discuss the 



essential points of the algorithm and its usage in our problem. 

Let us first recapitulate the definition of the discretization in the frequency uj and the 
imaginary time r; 

uoi = l-Auo, Ai = A{uJi)- ^uo, (/ = 1,2,3, ■■■,iV^), (3.26) 

= i-AT = i-a, Di = D{Ti), {i = 1,2,3,- ■ ■,N), (3.27) 

where Au is the pixel size of 0(10 MeV) in the frequency space (see Section 5.4 on how 
to choose Au;) and N^^ is the mesh size of O(IO^). Ar = a is the lattice spacing in the 
temporal direction. We have defined = Tmax/d — Tmin/o, + 1 ~ 0(10) as before. 
The extremum condition OQ/dAi = (eq.( |3.16| )) with eq.( |2.22D is written as 



^ Tmax/a. Qj^ 

-alog— = V Kii——, where Ku = K{Ti, ui). (3.28) 
nil ._ ^ , oDAi 

Since Ai is positive semi-definite, it is parametrized as 

Ai = m^expa^ (1 < / < A^"^). (3.29) 

Here a = {ai,a2,- ■ ■,a^^y is a general column vector in the A'^ = 0(10'^) dimensional 
space. However, a as a solution of ( p.28| ) turns out to be confined in the so called "singular 



subspace", whose dimension is Ns{< N <^ N^^). To show this, let us substitute ( p.29| ) 
into ( P^D , 



BL 

aa = K'-—, (3.30) 
oDa 



dL 

where is an x matrix and — is an A^ dimensional column vector. 

ODa _ 

The SVD of is defined by = IfEV^, where U is an A^^; x A^ matrix satisfying 
U^U = 1, V is an N X N matrix satisfying V^V = VV^ = 1, and S is an A^ x A'" 
diagonal matrix with positive semi-definite elements, C,i {i = 1,2, ...,N), which are called 
the singular values of ^j's may be ordered in such a way that > ^2 > ■ ■ ■ > 

^Ns > ^Af.+i = ■ ■ ■ = 0, where 

Ns = rank[ir*] < A^ < A^^. (3.31) 

The proof of the singular value decomposition is given in Appendix C for completeness. 
For the kernel such as {K*)h = exp(— u;/rj), the singular values are all non-zero but 
decrease exponentially as j increases. 
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The explicit form of SVD is written as 



UEV' 







V 



J 



\ 







.(3.32) 



Following Bryan , we define the Ns dimensional space spanned by the first Ng columns 
of U as the "singular space". The bases in this space are {ui,U2 - ■ ■ ,unJ\ with Ui = 
{uii,U2i, ■ ■ ■ yUN^iY ■ Then, one immediately observes from eqs.( |3.30| ) and ( |3.32D that a 
is in the singular space. This implies that a can be parametrized only by a set of A^^ 
parameters ■ ■ -bj^J as 



i=l 



Ui 



I.e., 



ai 



(3.33) 



i=l 



Therefore, owing to U^U = 1, eq. (|3.3CI|) reduces to 



ab = g = E' V" 



dL 



(3.34) 



where b is an Ng dimensional column vector, and S' and V are an Ng x A^^ matrix and 
an X A^s matrix obtained by restricting S and V to the singular space, respectively. 



In other words, they are defined by E'^j 



1,2, 



A^,) and V', 



1,2, 



iV, J = 1,2, 



,Ng) 



To solve (|3.34| ), the standard Newton method is used for each increment 6b as 
J 6b = —ab — o, with Jj,- = a/j,- + 

ob-i 



(3.35) 



where / is the identity matrix. By using the chain rule and the identity dA/db 
diag[A] U, ( p.35| ) is rewritten as 



[ (a + /i) / + MT] 6b = -ab - g, 



where 



M = E' 1/'* 



dD^dD^ 



V 



T = U" diag[A] t/'. 



(3.36) 



(3.37) 



with U' being defined by = U, 



l,2,---,iVs) 



Note that the Marquardt-Levenberg parameter /i is added to the diagonal element of 
the Jacobian matrix , so that the increment 6b at each iteration becomes small enough 
to guarantee the validity of the lowest order approximation used in the Newton method 
eq.( 3.35| ). As fi increases, the increment 6b generally decreases. The choice of fi is not 
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unique; here we follow ref. |^ and adjust so that the norm of 6 A defined by the metric 
gw = (l/Ai) 6ii' (see ( [A.9| ) in Appendix A) does not exceed the integrated default model: 

6 A diag[l/A;] 6A = 6b T 6b < c^, m , (3.38) 

1=1 

with c being a constant of 0(1). 

In our actual analysis, we use the SVD routine in ref. [0. Since the elements of the 
kernel vary many orders of magnitude, quadruple precision is necessary for obtaining 
reliable result of spectral functions at low frequencies. Then, at each iteration in ( |3.36| ), 
we start with fi = and increase /i by 10 multiples of 10~'^a until the norm condition 
with c = 0.2 is satisfied. Once this condition is fulfilled, the increment 6b is added to 
the temporary solution vector btemp, bnew = btemp + ^b. This process is iterated until 
\6Q/Q\ < 10"^ is achieved. 

Due to the correlation in the imaginary time direction in each Monte Carlo sample on 
the lattice, we must take into account the non-diagonal covariance matrix C defined by 
eq. (|3.9|) . In this case, the Bryan method is applied after the following transformations 

K K = R'^K, D = R-^D, (3.39) 

where i? is a matrix which diagonalizes the covariance matrix, R~^CR = diag [erf]. Since 
C is a real symmetric matrix, we take an orthogonal matrix for R, i.e., R~^ = R^. K 
and Z), instead of K and D, are taken as the kernel and the data, respectively. After the 
transformation, the likelihood function L is written as 

i^=2E \D^-Y:KaA,)| lal (3.40) 

Using this expression, we can carry out the Bryan method as in the case where the 
covariance matrix is diagonal. 

3.5 Uniqueness of the solution in MEM 



In order to show the uniqueness of the solution of (|3.16|) , we first prove the following 

proposition. 

Proposition: 

Consider a real and smooth function F with n real variables, namely F(xi, X2, ■ ■ ■ , x„) G R 
with {x\^X2i ■ ■ ■ ,a;„) G R". Suppose the matrix d"^ F / dxidxj is negative definite, i.e., 

E y^^-^Vj < (fo^ any G R - {0}), (3.41) 

then F has only one maximum if it exists. In other words, the solution of dF/dxi = 

(i = 1, 2, ■ ■ ■ , n) is unique if it exists. 

Proof: 

Assume that there are more than or equal to two solutions for dF/dxi = {i = 
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1,2, ■■■,n). Take any two of them xi and X2, and define an interpolation x{t) = 
xi + t{x2 — xi) = xi + ty, and G(t) = F{x(t)). From the assumption, dG{t)/dt is 
continuous and differentiable in [0, 1], and satisfies 



dG_ 



t=0 



dG 



0. 



(3.42) 



t=i 



Thus, from Rolle's theorem, there exists at least one t E (0, 1) such that 
d''G{t) A d^F 



dt^ 



yj 



0. 



(3.43) 



x=x(t) 



However, ( p.43|) contradicts (p.41|) . Hence there cannot be more than or equal to two 
solutions for dF/dxi = (i = 1, 2, ■ ■ ■ , n). If there is a solution for dF/dxi = {i = 
1, 2, ■ ■ ■ , n), it is the global maximum of F from eg. (|3.41| ). Thus the proposition is proved. 
(QED) 



We now proceed to prove that the solution of eq. (|3.16|) is unique and corresponds to 
the global maximum oi Q = aS — L if it exists. 

For an arbitrary A^^ dimensional non-zero real- vector z = {zi, Z2, - ■ ■ , zn^), aS satisfies 
that 



-a 



(3.44) 



where we have used < A; < oo and < a < oo. It is important to notice that the l.h.s. 
never becomes zero. 

On the other hand, from ( p.40| ). 



ii'=\ 



OAidAi 



^ z? 



-Zl' 



- ^ ^ < 0, with Zi = Y^ Kiizi. 



i=l 



a. 



(3.45) 



The l.h.s. of (|3.45 ) becomes zero in the direction where 5j = (i = 1, ■ ■ -A^). There are 
many such directions because the rank of K is at most A^, which is much smaller than 
A''^, (the dimension of the vector zi). This means that —L has a lot of flat directions, and 
there is no unique maximum of —L as a function of Ai. 

Once one adds ( |3.44| ) to ( |3.45| ), the solution of eq. (|3.16| ) becomes unique if it exists, 
because of the proposition just proved. The existence of aS (7^ 0) in Q = aS — L plays 
an essential role for this uniqueness. 



3.6 More on the covariance matrix 

The eigenvalue spectrum of the covariance matrix Cij is known to show a pathological 
behavior when Nconf is not large enough compared to A^. In fact, it is reported in ref. p2 
that the eigenvalue spectrum displays a sharp break when Nconf < N, i.e., some eigen- 
values are of similar magnitude, while the others are much smaller. This leads to the 
likelihood function L in (|3.40|) extremely large. Thus, it has been empirically preferred 
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to take Nconf > 2A^. In this subsection, we sliall show that tlie small eigenvalues found 
for Nconf < in ref. [^] are actually exact zeros. 



First, we rewrite the definition of the covariance matrix (|3.9|) as follows: 



N t N t 

C^J = EC:;= Y: dTdJ, (3.46) 

m=l m=l 

where is an x matrix defined for the m-th gauge configuration, and d^ = 

[NconfiNconf - l)]-'/HD"'{Ti) - D{t,)). 

Now, it is easy to show that Rank[Cj'^] = 1, since each column of reads 

df^i^d"^, d^, ■ ■ ■ , d^y and is thus proportional to the same vector [d"^Y. Let us consider the 
case when N^onf < N. Since Rank[A + B] < Rank[A] + Rank[i?] for arbitrary matrices, 
A and B, we obtain 

N r 
^ ' con J 

Rank[a,] < E Rank[C,7] = N^onf < N . (3.47) 

m=l 

As a result, the number of zero eigenvalues N^ero satisfies 

N.ero = N- Rank[a,] > A^ - N^onf . (3.48) 
Therefore, 

Nconf > N (3.49) 

is a necessary condition to avoid the pathological behavior of the eigenvalues of the co- 
variance matrix. 

As we shall see more in Section 5, a substantial number of sweeps for gauge configu- 
rations are inserted between measurements in the Monte Carlo simulation on the lattice. 
This is for carrying out measurements with minimally correlated gauge configurations. 
Thus, the direction of each is expected to be independent, and in most cases eq. (|3.49| ) 



is also a sufficient condition. In order to stabilize the behavior of the eigenvalues, how- 
ever, a condition such as Nconf > 2A^ may be imposed ||22[. In our lattice simulations 



discussed in Section 5, this condition is well-satisfied. We have also checked explicitly that 
our calculation is free of the pathological behavior and that the fiuctuations along the A^ 
principal axes obtained by diagonalizing the covariance matrix are well-approximated by 
Gaussian forms. 
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4 Analysis with mock data 



To check the feasibihty of the MEM procedure and to see the dependence of the MEM 
image on the quahty of the data, we made the following test using mock data. Key issues 
here are whether MEM can detect sharp peaks and fiat continuum simultaneously. To 
study this, we consider two cases: 

Schematic SPF: MEM using mock data obtained from a schematic spectral function 
with two Gaussian peaks; one is sharp and the other is broad. 

Realistic SPF: MEM using data obtained from a realistic spectral function parametrized 
to reproduce the e^e~ annihilation cross section into hadrons. 

In these tests, we have assumed T = and that the covariance matrix C is diagonal 
for simplicity. The non-diagonality of C, however, plays an important role in the case of 
the actual lattice QCD data. 



The basic strategy of the analysis is summarized as follows. 

(i) We start with an input image of the form Ai^ioj) = uj^pin{uj). uj^ is a factor expected 
from the dimension of the meson operators (ra = 2 in eq. (|2.15 )). Then, we calculate 



the mock data from eq. (|2.22|) at T = as 



D,n{n)= K{Ti,Uj) Ainiuj) du , (4.1) 

Jo 

where Umax may be chosen arbitrary large, but we simply take some reasonable 
number above which pin does not show appreciable variation. 

(ii) By taking Din{Ti) at discrete points and adding a Gaussian noise, we create a 
mock data Dmocki'^i)- To mimic the noise level of our lattice QCD data with the 
Dirichlet boundary condition in the temporal direction, the variance of the noise 
a(rj) is chosen as 

a(ri) = b X DiniTi) X ^ . (4.2) 

At 

Here the parameter b controls the noise level. and b are changed to see the 
sensitivity of the output image. 

(iii) We construct an output image Aout{0 < oj < ujmax) = ^"^Pouti.^) using MEM with 
the mock data Dmockiji)- Then, we compare Aout with Ai^. We need to introduce 
the resolution Aa; in the frequency space to perform the numerical analysis in MEM 
as explained in the previous section; Acj can be any small number. As a possible 
measure for the quality of the output image, we introduce a distance between pout 
and Pin as 

r = / [poutiuj) - pin{uj)] duj . (4.3) 
Jo 
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4.1 Schematic SPF 



To study how the maximum entropy method can reproduce simultaneously a sharp peak 
and a broad peak, we first consider the following schematic spectral function: 



Pi: 



1 



i=i,2 



v^(r,/2) 



,-(^-M,)V(r,/2)2 



(4.4) 



with (Mi,ri) = (1,0.01) GeV, (M2,r2) = (3,0.5) GeV. 

To mimic the lattice data discussed later, the parameters for the analysis here are 
chosen to be uJmax = 6 GeV, Au = 10 MeV, Tmin = Ar = 0.085 fm. To see how Aout is 
improved as the quality of data increases, = Tmax/Ar and b are changed within the 
intervals, 10 < iV < 30 and 0.0001 < b < 0.01. As for m, we take the form m{uj) = itiquo'^, 
where mg is chosen to be 0.36 so that pin{uj)uj'^duj = niQU^duj is satisfied. 

In Fig.|I|, a comparison of Poutil^) (the solid line) and Pinij^) (the dashed line) is shown 
for various combinations of and b. The distance r defined in (|4.3|) is also shown in the 
figure. Increasing A^ and reducing the noise level b lead to better SPFs closer to the input 
SPF as is evident from the figure. MEM reproduces not only the broad peak but also 
the sharp peak without difficulty. This is because the entropy density in (|3.12| ) is a local 
function of uj without any derivatives and hence does not introduce artificial smearing 
effect for SPF in the cj-space. Fig.|I| also shows that, within the range of parameters 
varied, decreasing b is more important than increasing A^ to obtain a better image. 
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Figure 1: Input SPF with two Gaussian peaks (the dashed lines) and output SPF obtained 
by MEM (the solid lines) for different values of data-points A^ and noise level b. 

In Fig.^ the probability distribution P[a\DHm] (with an approximation discussed 
in Step 2 in Section 3.3) used to obtain the final image is shown for three different 
combinations of A^ and b. The distribution tends to become peaked as A^ increases and b 
decreases. 
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Figure 2: Probability distribution P[a\DHm] for three different sets of (A^, b) in the case 
of Schematic SPF. 

4.2 Realistic SPF 

As an example of realistic spectral functions, we study SPF in the charged p-meson 
channel. The isospin symmetry relates it to the e~^e~ annihilation data into hadrons in 
the isospin 1 channel. We take a relativistic Breit-Wigner form as a parametrization of 
SPF in this channel 0: 



f ^ 2 

IT 



p2 r^mp ,J_f-i,^^ 1 



(4.5) 



P (cj2 - m2)2 + r2m2 Svr V tt / 1 + 6(^0-"^)/^ 

The pole residue Fp is defined: 

(0|J7^m|p) = V2Fpmpe^ = v^/pm^e^, (4.6) 

where is the polarization vector. In the vector dominance model, the dimensionless 
residue fp is related to the pirn coupling gp^.^^ as fp = l/gpTr^- 

To get the correct threshold behavior due to the p — tttt decay, we take the following 
energy- dependent width: 



( 4m2 \ 

r,(.;) = ^ m J 1 - e{^-2m;). (4.7) 



The empirical values of the parameters are 

mp = 0.77 GeV, = 0.14 GeV, 
gp^^ = 5.45, (4.8) 
ooo = 1.3 GeV, 6 = 0.2 GeV, 

where we have assumed the vector dominance and ag is taken to be 0.3 independent of cu 
for simplicity. 

As for m, we take the form m{uj) = ■rriQUj'^ motivated by the asymptotic behavior 
of Pin{oj) in ( [4.5[ ). mo is chosen to be 0.0257, which is slightly smaller than 0.0277 ex- 
pected from the large uj limit of ( [4.5|) , Pm(^ — *■ oo) = (2/7r)(l/87r)(l + a^/vr). The 
same parameters as those for the schematic SPF are chosen: uJmax = 6 GeV, Aa; = 10 
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MeV, r„ 



At = 0.085 fm. = Tmax/Ar and h are changed within the intervals, 



10 < < 30 and 0.0001 < 6 < 0.01. 

In Fig.^ the mock data Dmock{Ti) obtained from ( ^.1| ) and ( [4. 51 ) are shown for the 
noise parameter b = 0.001. The hnear slope of logDmock{Ti) for large Tj is dictated by the 
lowest resonance, while the deviation from the linear slope for small Xj originates from the 
continuum in SPF. 
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Figure 3: Mock data Dmock{Ti) in the case of the realistic SPF with Gaussian error 
attached. 
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Figure 4: Input SPF with a resonance + continuum (the dashed lines) and output SPF 
obtained by MEM (the solid lines) for different values of N and b. 

In Fig.§, a comparison of Poutil^) (the solid line) and Pinif^) (the dashed line) is shown 
for various combinations of N and b. The distance r defined in ( |4.3| ) is also shown in the 
figure. Increasing and reducing the noise level b lead to better SPFs closer to the input 
SPF as is evident from the figure as in the case of the schematic SPF. 
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In Fig.^ the probability distribution P[a\DHm\ (with an approximation discussed 
in Step 2 in Section 3.3) used to obtain the final image is shown for three different 
combinations of N and h. The narrower distribution is obtained when the data quality is 
better. 



0.0025 




a (GeV"') 

Figure 5: Probabihty distribution P[a\DHm\ for three different sets of {N, h) in the case 
of the realistic SPF. 




Figure 6: Default model dependence of the output image Pout{^) (solid lines) for = 25 
and b = 0.001. The input images and the default models are shown by the dashed 
and dotted lines, respectively. Error bars for three different frequency regions are also 
attached. The middle figure is the same as that in Fig.^. 
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Finally, to see the sensitivity of the results against the change of mo, pouti,^) is shown 
in FigJ^ for five different values of the default model together with the error bars. = 25 
and h = 0.001 are chosen. Short dashed and long dashed lines correspond to the default 
model and the input SPF, respectively. The solid lines show the output image after the 
MEM analysis, ttiq = 0.257 is chosen in the middle figure which is the same with that 
in Fig.^ while the other four figures correspond to the default models for 0.25mo, 0.5mo, 
2mo and 4mo. Even under the factor 4 variation of the default model, the resultant SPFs 
show the peak + continuum structure. However, as the default model deviates from the 
expected asymptotic value, the SPF starts to have a "ringing" behavior. 

Now, the error analysis discussed in Step 3 in Section 3.3 can tell us whether the 
ringing structure seen, e.g., in the 4mo case (the upper right figure) corresponds to a real 
resonance or is the artifact of the maximum entropy method. The horizontal position and 
length of the bars in Fig.|^ indicate the frequency region over which the SPF is averaged, 
while the vertical height of the bars denotes the uncertainty in the averaged value of the 
SPF in the interval. Implications of the error bars in Fig.^ are twofold; (i) the ringing 
images for Aitlq and 0.25mo are statistically not significant, and (ii) the combination of 
the best SPF and the best default model may be selected by estimating the error bars 
(Step 4 in Section 3.3). In the present case, for the given data, the middle figure should 
be chosen to be the best one. 
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5 Analysis with lattice QCD data 



5.1 Lattice basics 

In lattice gauge theories the Euchdean space-time is discretized into cells. In our 
simulations, we use an isotropic hypercubic grid, namely the lattice spacing a is the same 
in all directions. In the following, we assume that the gauge group is SU(3) with QCD in 
mind. The fermion field il){x) is defined at grid sites, while the gluon field is defined on 
links connecting adjacent pairs of grid sites (Fig.|^). On the links connecting x to x ± /i 
= a), we define ?7±^( clS 



U^{x) = exp [iagAf,{x)] , f/„ 



exp [—iagA^{x 



(5.1) 



f/^(x) is an element of the SU(3) group, g is the gauge coupling constant, and is the 
gauge field. The local gauge transformation for ip^x) and U^{x) on the lattice reads 



■ip{x) 
U^{x) 



V{x)il>{x), il>{x) il>{x)V\x) 
V{x + fi)U^{x)V\x), 



(5.2) 



where V{x) is an SU(3) matrix defined at each site x. 
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Figure 7: Lattice variables. Quarks (gluons) are defined at the sites (on the links). 



Using the plaquette variable defined by 

W{x,fi, 0) = U^u{x + 0)U-^{x + fi + u)Uu{x + ij)U^{x), 

the single plaquette gluon action is written as 

X 6 
Sp = (3iat^^^eTi-{l-W{x,fi,i))) , with (3iat = —■ 

X ii<u 9 

In the naive continuum limit a — 0, Sp approaches the continuum gauge action 
1 



Scout = ^ I d'xiF^^ix)y. 



(5.3) 



(5.4) 



(5.5) 
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For fermions, we use the Wilson quark action defined by 



x) 



x,y 



(l + 7/.)^^(y) 



= J2^{x)D^ix,y;U)^lj{y). (5.6) 

Here the flavor indices for fermions are suppressed. 7^ satisfies {7/x,7i/} = 25^j,. is 
taken for both positive and negative directions of /i together with a convention 7_^ = — 7/^- 
K denotes the hopping parameter, while the quark mass rrig is defined as 

11 , , 

for a fixed jSiat- Here KciPiat) is the critical hopping parameter, at which the pion becomes 
massless. 

The action for the whole system Siat is given by the sum of 5*^ and S^, 

Slat = Sp + Sw (5.8) 

The expectation value of a physical observable 0[U,ip, ^] is given in terms of path integrals 
over the link variables U^{x) and fermion fields iplx) and ip{x), 

{0[U, = 1 /" [dUd^dtp] e-^'"* 0[U, tp, (5.9) 

Zj J 

where Z is given by 



Z 



[dUdipdtp]e--'^"^' = J [dU]e--'^^ detD^{x,y; U). (5.10) 



In the quenched approximation, detDw{x,y]U) is set to 1. Physically, this corresponds 
to ignoring the effect of virtual quark loops. This approximation simplifies numerical 
simulations by factor 1000 or so, since it is very time consuming to calculate the deter- 
minant of a huge matrix such as D^{x,y; U). In the quenched simulations, the first step 
is to generate an ensemble of gauge configurations U^{x) with the weight e""^*". A typical 
method for that purpose in SU(3) gauge theory is the pseudo heat-bath method combined 
with the over-relaxation method (for more details, see, e.g., p9|). 

5.2 Lattice parameters 

We take the open MILC code with minor modifications and perform simulations with 
the single plaquette gluon action + Wilson quark action in the quenched approximation 
on a Hitachi SR2201 parallel computer at Japan Atomic Energy Research Institute. The 
basic lattice parameters in our simulations are 

6 

couplmg strength fiiat = ^ = 6.0, 

lattice size Nl x = 20^ x 24, 
hopping parameter n = 0.153, 0.1545 and 0.1557. (5-11) 
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The Dirichlet boundary condition (DB) in the temporal direction, which is defined by 
Unix) = for the last temporal links, is employed to have as many temporal data points as 
possible. In the spatial directions, the periodic (anti-periodic) boundary condition is used 
for the gluon (quark) as usual. Gauge configurations are generated by the pseudo heat- 
bath and over-relaxation algorithms with a ratio 1 : 4. Each configuration is separated 
by 200 sweeps. The number of gauge configurations used in our analysis is Nconf = 161. 

We have also carried out simulations with the periodic (anti-periodic) boundary condi- 
tion for the gluon (quark) in the temporal direction on the 20'^ x 24 lattice with f3iat = 6.0, 
and on the 40^ x 30 lattice with Piat = 6.47 on CP-PACS at Univ. of Tsukuba. Detailed 
MEM study in those cases will be reported elsewhere El . 



To calculate the two-point correlation functions, we adopt a point-source at (r, x) = 
(0, 0), and a point-sink at time r with the spatial points averaged over the spatial lattice 
to extract physical states with vanishing three- momentum (see eq. ( p.22| )). The following 
local and fiavor non-singlet operators are adopted for the simulations: 

scalar (S) : du, pseudo-scalar (PS) : d'y^u, (5-12) 

vector (V) : d'jf^u, axial-vector (AV) : ^7^75^. (5.13) 

In the V and AV channels, the spin average is taken over the directions fi = 1,2,3 to 
increase statistics. Simulations for spin 1/2 and 3/2 baryons have been also carried out. 
Since special considerations are necessary for the decomposition of the SPFs in the baryon 
channels, we shall report the results in a separate publication |^ . 



Note here that the use of the point source and sink is essential for obtaining good 
signals for the resonance and continuum in the spectral function simultaneously. First of 
all, the SPF defined with the point source and sink for the vector channel is directly related 
to the experimental observables such as the i?-ratio and the dilepton production rate as 
discussed in Section 2. Besides, there are two- fold disadvantages to use smeared sources 
and sinks in the MEM analysis: First of all, it is difficult to write down a simple sum rule 
such as eq. ( p.22|) for the smeared operators. Secondly, the coupling to the excited hadrons 



becomes small for such operators and one loses information on the higher resonances and 
continuum. 

For the temporal data points to be used in the MEM analysis, we take T^in/a = 1 
and Tmax/d = 12. The latter is chosen to suppress the error from the Dirichlet boundary 
condition. In fact, we found that the statistical error of our data is well parametrized 
by formula ( [4. 21 ) with a slope parameter 6 up to r/a = 12 and that the error starts to 
increase more rapidly above r/a = 12. 

5.3 Spectral functions and their asymptotic forms 

We introduce the dimensionless SPFs, p{uj), as follows: 

M=l,2,3 

Ps psvAvi^) ^'^^ defined in such a way that they approach finite constants when u ^ 00 
in the continuum limit (a — > 0) as predicted by perturbative QCD (see eq.( |2.15[ )). 
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Because of ( p.8[ ) and (|2.9|), we have 



Ps,ps,v,Avi^) = -Ps,ps,v,Avi-^) > 0, with ^ > 0. (5.15) 

The SPFs on the lattice in the chiral hmit should have the following asymptotic form 
when uj is large and a is small; 

» 1 GeV) ^ (l + [^^)\ (5,16) 

where j = S, PS, Vand AV. The first two factors on the r.h.s. are the qq + qqg con- 
tinuum expected from perturbative QCD. The third factor contains the non-perturbative 
renormalization constant, Zj{pa), of the lattice composite operator. The renormalization 
point /i should be chosen to be the typical scale of the system such as 1/a. 

Tij and calculated perturbatively in the continuum QCD are shown in Table 



In [Q, Zj{fia = 1) in the chiral limit has been calculated numerically on the lattice 
with Plat = 6.0 for PB (the periodic (anti-periodic) boundary condition in the temporal 
direction for the gluon (quark)), which is summarized in Table |l]. Zj for DB and those 
for PB are in principle different. This is because, in our simulation, the hadronic source 
is always on a time slice where DB is imposed and Zj detects the the boundary effect. 
Our measurement of the decay constant fp confirms this as shown later in section 5.5.1. 
Therefore, we use Zj listed in Table [l| only to guide our default model m. 



J 




r2j 


Zj{fia = 1) 


Pj{uj=fx=l/a) 


s 


3/2 


11/3 


0.77 


0.80 


PS 


3/2 


11/3 


0.49 


2.00 


V 


1 


1 


0.68 


0.59 


AV 


1 


1 


0.78 


0.45 



Table 1: Constants for spectral functions in the asymptotic region, rij and r2j are taken 
from Zj{fia = 1) for j3iat = 6.0 in the chiral limit Hc = 0.1572 is taken from [Q, 



where PB is used in the temporal direction. For the estimates in the last column, we use 
a,(/i ~ 2 GeV) = 0.3 |5|, = 0.1572 and ~ 2 GeV for piat = 6.0. 



So far we have assumed that the spatial momentum k vanishes due to the spatial 
integration on the lattice. However, it is not exactly the case for any finite set of gauge 
configurations. In most cases, the finite k error is harmless as far as k is small enough. 
However, it is potentially dangerous in the AV channel, where the SPF for finite = 
(a;, k) is written as 

A^^{LJ,k) = {k^K - k'^g^y) p^{uj,k) + k^k^ p^{uj,k). (5.17) 

Here Pj, is the transverse spectral function, which does not have contamination from 
the pion pole. When k = 0, reduces to p^y{uj > 0) > defined in ( |5.14| ). On 



the other hand, is the longitudinal spectral function, which contains the pion pole: 
p^(a;,k) = 2fl5{k'^ — m^) + ■ ■ ■ with being the pion decay constant. Then, if small 
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amount of k remains, there is possible contamination from tlie low-mass pion pole to the 
spin-averaged SPF: 

\ V(u;,k^O)^o.X,M + kV,M. (5.18) 

At=l>2,3 

In the actual lattice data, this effect appears in D{t) for large r. A possible way to 
subtract the contamination of is to carry out the measurement with finite k. 



5.4 Discretization in uj space 

In the MEM analysis, we need to discretize the cu-space into pixels of an equal size Au 
and carry out the integration, ( p.22| ), approximately. The upper limit of Au is determined 
by the requirement that the discretization error of the kernel K in ( p^.22| ) is small enough, 
namely, r ■ Au <^ 1, which reads 

Alj < < -. (5.19) 

Thus r~^^ = (A^^a)"^ = (24a) ~ 90MeV is the upper bound of Aoj on our lattice. 

Also, to obtain the good resolution of SPF, one needs to have small enough Auo so 
that — yli|/|y4i_|_i + Ai\ ^ 1. In the maximum entropy method with singular value 
decomposition, there is no problem in choosing an arbitrary small value for Auj. In fact, 
as we have explained, the maximum search of Q is always limited in the Ns dimensional 
singular space, and Ng is independent of N^. Therefore, increasing N^^ or decreasing Auj 
does not cause any numerical difficulty. In the following, we adopt Aa; = 9.5 MeV, which 
satisfies ( p.l9[ ) and simultaneously gives a good resolution to detect sharp peaks in the 



spectral function. 

Aside from Auj, we also need to choose the upper limit for the uo integration, uJmax- 
Since this quantity should be comparable to or larger than the maximum available mo- 
mentum on the lattice, we choose Umax = 7.1 GeV > n/a 6.9 GeV. We have checked 
that larger values of uJmax do not change the result of ^(0;) substantially, while smaller 
values of uJmax distort the high energy end of the spectrum. The dimension of the im- 
age to be reconstructed is N^^ = uJmax/Auj ~ 750, which is in fact much larger than the 
maximum number of data points, A^ = 24, available on our lattice. 

5.5 Results of MEM analysis 
5.5.1 Pseudo-scalar and vector channels 

Let us first consider the PS and V channels. The lattice data D{Ti) in these channels are 
shown in Fig.|^. For the MEM analysis, we take the formula ( ^.22| ) with (3 = 1/T = 00, 
since we use DB (the Dirichlet boundary condition). Also, Tmin/O' = 1 has been taken to 
utilize the information available on the lattice as much as possible, while Tmax/O' = 12 has 
been chosen to suppress errors caused by DB as discussed in the previous section. 

Shown in Figs.P(a) and ^(b) are the reconstructed SPFs in these channels for different 
values of k. Umax IS taken to be 7.1 GeV and Au = 9.5 MeV. We have used m = rriouj^ 
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with rriQ = 2.0 (0.6) for the PS (V) channel motivated by the perturbative estimate 
mo ~ p{uj ^ 1 GeV) in eq. ( ^.16| ) with Table |I[ The sensitivity of the reconstructed SPFs 
on the variation of mo will be discussed later. 

Spectral functions thus obtained in the PS and V channels shown in Figs.|^ and 
have a common structure: a sharp peak at low-energy, a less pronounced second peak, 
and a broad bump at high-energy. Let us discuss those structures in detail below. 
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Figure 8: Lattice data -D(rj) in the PS and V channels with statistical errors. 
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Figure 9: Reconstructed image Poutil^) for the PS (a) and V (b) channels. The solid, 
dashed and dash-dotted lines are for /t = 0.1557, 0.1545 and 0.153, respectively. For the 
PS (V) channel, mo is taken to be 2.0 (0.60). Umax is 7.1 GeV. 



(i) To make sure that the lowest peak for each k in the PS (V) channel corresponds to 
the pion (the p- meson) , we extract m^r and irip in the chiral limit with the following 
procedure. First we fit the lowest peak by a Gaussian form and extract the position 
of the peak. Then, the linear chiral extrapolation is made for (m^ro)^ and iripa. 
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The results together with a using the input rup = 0.77 GeV are summarized in 
the column "MEM continuum kernel" in Table ^. They are consistent with those 
determined from the standard analysis using the asymptotic behavior of -D(Tj) in 
the interval 10 < Tj/a < 15, which are shown in Table ^ in the column "asymptotic 
analysis" . Also, our results are consistent with those from the asymptotic behavior 
of D{t) obtained by the QCDPAX Collaboration on a 24^ x 54 lattice with Piat = 6.0 
H. They give = 0.1571, m^a = 0.331(22) and a'^ = 2.33(15) GeV. 

Although it is quite certain that the lowest peaks for each k in Figs.|^(a) and |^(b) 
correspond to vr and p, the widths of these peaks in the quenched approximation do 
not have physical meaning, since they are artifact caused by the incompleteness of 
the information contained in lattice data. Nevertheless, the integrated strength of 
the peak corresponds to the physical decay constant of the mesons. For example, the 
dimensionless decay constant fp defined in eq.( [1.6| ) is related to the lattice matrix 
element as 

{0\{d^,uU\pi^ = 0)) = v^e^/pm^-J— . (5.20) 

AKZjy 



Therefore, it is further related to the spectral integral near the resonance as 

2 

/pole' - ' ' VZkZv.' 



j pA.W = 2fl,nl [^)\ (5,21) 



where the suffix ^^pole" implies the integral over the lowest isolated peak. This 
integral can be carried out numerically for each k, and the linear extrapolation to 
the chiral limit has been done. The result for fp/Zy is given in the 2nd row of 
Table ^ together with that obtained from the asymptotic analysis. The agreement 
is satisfactory. As we have mentioned before, Zy in the DB case is not necessary 
equal to that in the PB case. Therefore, we cannot predict fp from our data alone. 
Instead, we extract Zy in the DB case by comparing our measurement of fp/Zy with 
experimental value of fp. This leads to Z^^ ~ 0.38 < Z^^ ~ 0.68. If we can have 
larger lattice and can place the hadronic source and sink far from the boundary, the 
difference between Z^^ and Z^^ should become small. 

As for the second peaks in the PS and V channels, the results of the error analysis. 



which are shown in Fig.[Iy, indicate that their spectral "shape" does not have much 
statistical significance, although the existence of the non-vanishing spectral strength 
is significant. With this reservation, we fit the position of the second peaks and make 
linear extrapolation to the chiral limit. The results are summarized in the 3rd row 
of Table § together with the experimental values. Since there exist two excited-p 
experimentally observed, p(1450) and p(1700), we quote both values for rrip'/mp in 
Table |. 

One should note here that, in the standard two-mass fit of D{t), the mass of the 
second resonance is highly sensitive to the lower limit of the fitting range, e.g., 
m2"'^/mp = 2.21(27) (1.58(26)) for r™„/a = 8 (9) in the V channel with (3iat = 6.0 



. This is because the contamination from the short distance contributions, which 
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dominates the correlators at r < Tmm, is not under control in such an approach. 
On the other hand, MEM does not suffer from this difficulty and can utilize the 
full information down to Tmm/a = 1- Therefore, MEM opens up a possibility of 
systematic study of higher resonances with lattice QCD. 

(iii) As for the third bumps in Fig.^, the spectral "shape" is statistically not significant 



as is discussed in conjunction with Fig. 10. They should rather be considered to 



be a part of the perturbative continuum instead of a single resonance. Fig.| and 
Fig.|l^ given later show that SPF decreases substantially above 6 GeV, namely MEM 
automatically detects the existence of the momentum cutoff on the lattice ~ vr/a. 
It is expected that MEM with the data on finer lattices leads to larger ultraviolet 
cut-offs in the spectra. Our preliminary analysis on a finer lattice (40^ x 30 lattice 
with Plat = 6.47) in fact indicates that this is the case HTI . 





asymptotic analysis 
(10 < Til a < 15) 


MEM 
continuum kernel 

(1 < Ti/a < 12) 


MEM 
lattice kernel 

(1 < n/a < 12) 


Experiments 


rUpa 
(GeV) 


0.1572(1) 
0.343(10) 
2.24(7) 


0.1570(3) 
0.348(15) 
2.21(10) 


0.1569(1) 
0.348(27) 
2.21(17) 




fp/Zv 


0.48(3) 


0.520(6) 




/p =0.198(5) 


rripi/mp 




1.88(8) 
2.44(11) 


1.74(8) 
2.25(10) 


1.69(13) 
1.90(3) or 2.21(3) 



Table 2: A comparison of the MEM analysis and the asymptotic analysis of the same 
lattice data obtained on the 20'^ x 24 lattice with (3iat = 6.0. The experimental value of 
fp obtained from the decay p e^e~ and the masses of p' and tt' are also shown [^] . 
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Figure 10: Pout{'^) in the V and PS channels for k = 0.1557 with error attached. uJmax = 7.1 
GeV. 
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In Fig.|TO|, we show the MEM images in the V and PS channels for k = 0.1557 with 
errors obtained in the procedure discussed in Step 3 in Section 3.3. The meanings of the 
error bars are the same as those in Fig.|^. Namely, the horizontal position and length of 
the bars indicate the frequency region over which the SPF is averaged, while the vertical 
height of the bars denotes the uncertainty in the averaged value of the SPF in the interval. 



The small error for the lowest peak in Fig.|TO| supports our identification of the peak 
with p. Although the existence of the non-vanishing spectral strength of the 2nd peak 
and 3rd bump is statistically significant, their spectral "shape" is either marginal or 
insignificant. Lattice data with better quality are called for to obtain better SPFs. 




Figure 11: P[a\DHm] for the PS and V channels. 



In Fig.pA], P[a\DHm\ (with the approximation discussed in Step 2 in Section 3.3) is 
shown for the PS and V channels in the case of k = 0.1557. We have found that the SPFs 
obtained after averaging over a and those at a = a have negligible difference although 
the distribution of a spreads over the range ^ 1 in Fig.|TT[ 

The sensitivity of the results under the variation of uimax is shown in Fig.lT2|, where 



^max = 8.5 GeV is adopted. Comparison of Fig.|T2|and Fig.|| shows no appreciable change 
of SPFs under the variation of uJmax as long as ujmax > vr/a. 

We have also checked that the result is not sensitive, within the statistical significance 
of the image, to the variation of itlq by factor 5 as shown in Fig.|T3|. The default model 
dependence is relatively weak compared to the case of the mock data shown in Fig.^ 
partly because the off-diagonal components of the covariance matrix for the lattice data 
are not negligible and stabilize the final images. 



^It is in order here to make some comments on the difference in the figures shown here and those in 
our previous publications ||l6| . The lattice data used for the MEM analyses are exactly the same in both 
cases. In a^^ = 2.33 GeV (which is obtained on 24^^ x 54 lattice in the first reference of [Q) was 
used to set the scale, while in this article we use — 2.21 GeV obtained from the p-meson mass in 
our MEM analysis of the data on the 20'^ x 24 lattice. This leads to a simple rescaling of the factor 0.95 
for dimensionful quantities such as LOmax and Auj from those in [^6|. Also, we mentioned in []l6| that 
the averaged heights of the high energy continuum of SPF in the V and PS channels are consistent with 
the perturbative prediction in Table |l|. This statement is misleading, since Zj"^ (which are shown in the 
table) are different from Zj^^ to be used in our analysis. 
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Figure 12: SPFs in the V and PS channels for Umax = 8.5 GeV. The sohd, dashed and 
dash-dotted hnes are for k = 0.1557, 0.1545 and 0.153, respectively. 
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Figure 13: The default model dependence of the output image in the vector channel for 
^max = 7.1 GeV. The solid line corresponds to mo = 0.6 (the same value with that taken 
in Fig.H and Fig.^, while the long (short) dashed line corresponds to 5mo (mo/5). 



5.5.2 Scalar and axial- vector channels 

In Fig.Il, the spectral functions for the S and AV channels are shown. In those channels, 
we have used Tmax/O' = 12 (10) for the S (AV) channel. Smaller Tmax is chosen in the AV 
channel to avoid the pion contamination that contributes to the correlator at large r as 
discussed in Section 5.3. Although the peak + continuum structure is seen in Fig.|l^, SPF 
does not change in a regular way under the variation of the quark mass, which prevents us 
from making chiral extrapolation of the peaks. More statistics and also better technique 
are needed to obtain SPFs with good quality in these channels. 
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Figure 14: SPF for S and AV channels. The sohd, dashed and dash-dotted hnes are for 
K = 0.1557, 0.1545 and 0.153, respectively. 



5.5.3 Lattice versus continuum kernel 

In the MEM analysis presented so far, we have relied on the spectral representation derived 
in the continuum limit, 



D{t) = / K{t,uj) A{iu) duj. (5.22) 

JO 

Here the kernel at T = is related to the free boson propagator with a mass u: 



e 

K{t, u) = = 2uj I ^. (5.23) 

Now, to study the systematic error from the finite lattice spacing in extracting SPF, 
one may artificially define a "lattice spectral representation" as 

D(r)= K^''\t,uj) A^^Uuj) du, (5.24) 

JO 

where K'-"-^ is a "lattice kernel" at T = defined through the lattice boson propagator 
with a mass uj: 

/.-'(r..;a)^2./ . (5.25) 
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j^iat ^^^^ j^iat approach K and A, respectively, in the continuum limit. Therefore, by 
taking the same lattice data D{t) on the l.h.s. of ( |5.22|) and ( |5.24|) and by comparing the 
resultant A{uj) and A'-"'^{uj) in the MEM analysis, one can estimate a part of the systematic 
error from the finiteness of the lattice spacing. 

The difference of iiT'"* and K becomes substantial for large tu, which is shown in 
Fig.[l^. Therefore, the finite a error appears typically at large u. The spectral functions 
in the V and PS channels obtained with the lattice kernel are shown in Fig.[T^. In Table ^ 
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the masses extrapolated to the chiral hmit are shown in the column "MEM lattice kernel" . 
{fp is not shown in the table, since it is difficult to isolate the p-pole unambiguously from 
Fig.[T^.) The results are consistent with those of "MEM continuum kernel" within error 



bars. In Fig.|T^, the SPFs with error attached are shown for the PS and V channels. 
Taking into account those errors, it is hard to distinguish the images obtained by K and 
fC'"* with the present lattice data, although the SPFs in Fig.|l^ are flat compared with 
those in Fig.|To|. 
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Figure 15: Comparison of the lattice kernel fC''^* and the continuum kernel K near the 
highest frequency on the lattice lj = 3 /a. 
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Figure 16: SPFs for the V and PS channels obtained with the lattice kernel. The solid, 
dashed and dash-dotted lines are for k = 0.1557, 0.1545 and 0.153, respectively. 
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Figure 17: SPFs for the V and PS channels obtained with the lattice kernel with error 
attached for k = 0.1557. 
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6 Summary and concluding remarks 



In this article, we have examined the spectral function (SPF) in QCD and its derivation 
from the lattice QCD data. The maximum entropy method (MEM), which allows us to 
study SPFs without making a priori assumptions on the spectral shape, turns out to be 
quite successful and promising even for limited number of lattice QCD data with noise. In 
particular, the uniqueness of the solution and the quantitative error analysis make MEM 
superior to any other approaches adopted previously for studying SPFs on the lattice. 
Also, the singular value decomposition (SVD) of the kernel of the spectral sum rule leads 
to a very efficient algorithm for obtaining SPFs numerically. 

By analyzing mock data, we have tested that the MEM image approaches the exact 
one as the number of temporal data points is increased and the statistical error of the 
data is reduced. Even with the lattice QCD data at T = obtained on the 20^ x 24 lattice 
with the lattice spacing a ^ 0.09 fm, MEM was able to produce resonance peaks with 
correct masses and the continuum structure. The statistical significance of the obtained 
spectral functions has been also analyzed. Better data with smaller a and larger lattice 
volume will be helpful for obtaining SPFs with smaller errors. 

MEM introduces a new way of extracting physical information as much as possible 
from the lattice QCD data. Before ending this article, we list open problems for future 
studies. Some of them are straightforward and some of them require more work. 



Technical issues: 



The number of temporal points and the lattice spacing a are the crucial quantities 
for the MEM analysis to be successful. In this article, we have fixed = 12 and 
a = 0.09 fm for extracting SPF from the lattice data. However, it is absolutely 
necessary to study A^ and a dependence of the resultant SPFs. For this purpose, 
we have collected 160 gauge configurations on CP-PACS at Univ. of Tsukuba with 
40^ X 30 lattice and Aat = 6.47 (a ~ 0.05 fm) with the periodic boundary condition. 
The detailed MEM analysis of the data will be reported elsewhere . 



(2) We have made chiral extrapolation only for the peaks of SPFs but not for the 
whole spectral structure in this article: In fact we have found that neither the 
direct extrapolation of the MEM image Aouti^^) nor the extrapolation of D{t) and 
C works in a straightforward manner. This is an open problem for future study. 

(3) Although we have shown that one can select an optimal default model by vary- 
ing m(cij) and estimating the errors of SPFs, the variation was within an assumed 
functional form such as m{uj) = mQUj"' motivated by the asymptotic behavior of 
the spectral functions in QCD. How to optimize the default model given data in a 
systematic way should be studied further. (For a possible procedure by estimating 
P[m\DH], see ref.|22|.) 
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Physics issues: 



(4) Further studies on the scalar and axial-vector channels are necessary to explore the 
applicability of MEM for not-so-clean lattice-data. 

(5) MEM analysis for baryons is quite useful in extracting information on excited 



baryons and their chiral structure 



(6) Applications to heavy resonances such as the glueballs, and charmed/bottomed 
hadrons will be an ideal place for MEM, since one can extract ground state peaks 
with limited number of data points obtained at relatively short distances. For the 



study of charmonium on 40 x 30 lattice with Piat = 6.47, see 41 



(7) When one studies a two-point correlation function of operators Oi and O2 with 
Oi ^ O2, one encounters SPF which is not necessarily positive semi-definite. Typical 
examples are the meson and baryon mixings such as p-u, tt^-t], rj-rj', glueball-gg and 
A-E" (see e.g. P^). The generalization of the Shannon- Jaynes entropy to non 



positive-definite images is possible [48|, and it is an interesting future problem to 



study hadron mixings from the generalized MEM. 

(8) Full QCD simulations combined with MEM may open up a possibility of first prin- 
ciple determination of resonance widths such as p(770) — 2it, (t(400 — 1200) 2it 
and ai(1260) Sir. Also it serves for unraveling the structure of the mysterious 
scalar meson "cr" [^. 

(9) The long-standing problem of in-medium spectral functions of vector mesons (p, u, 
0, J/'ip, T, ■ ■ -, etc ) and scalar /pseudo-scalar mesons (a, vr, ■ ■ -, etc.) can be studied 
using MEM combined with finite T lattice simulations. The in-medium behavior of 
the light vector mesons |10|, ^ |2^, ^ and scalar mesons [0, ^ is intimately 



related to the chiral restoration in hot /dense matter, while that of the heavy vector 
mesons is related to deconfinement |[Tl|, |5^, Q . An anisotropic lattice is necessary 
for this purpose to have enough data points in the temporal direction at finite T 
^5[ . (See also a recent attempt on the basis of NRQCD simulation 0].) Also, it 



is interesting to study SPFs with finite three-momentum k, since u and k enter in 
the in-medium SPFs as independent variables. 

(10) Possible existence of non-perturbative collective modes above the critical temper- 



ature of the QCD phase transition speculated in |5^, |58[ may also be studied 
efficiently by using MEM, since the method does not require any specific ansatze for 
the spectral shape. Also, correlations in the diquark channels in the vacuum and in 
medium are interesting to be explored in relation to the q-q correlation in baryon 
spectroscopy and to color superconductivity at high baryon density . 
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A Monkey argument for entropy and prior probabil- 
ity 

Throughout Appendix A and B, we shall use the common notation f{x) for the image 
instead of A{u) used in the text. 

What we need in the maximum entropy method is the prior probability P(/ G V^), 
namely, the probability that the image f{x) is in a certain domain V. It can be generally 
written as 

Pif ^V) = [W] ^aS{f)), (A.l) 



Zs{a) Jv 

where a is an arbitrary constant defined for later use, and Zs{(y) is a normalization factor. 
We assume that $ is a monotonic function of the would-be entropy S{f), so that the most 
probable image / is obtained as a stationary point of S{f). 

The so-called "monkey argument" , which is based on the law of large numbers, can 



determine the explicit form of $ and S{f) (see, e.g., [^T], ^ |3^, Q). Here we shall 
recapitulate the essential part of this derivation. 

Let us first discretize the basis space x into N cells. Correspondingly, f{x) is dis- 
cretized as (1 < i < N). Suppose a monkey throws M balls. M is assumed to be large. 
Define rii as the actual number of balls which the i-th cell received. Also, the i-th cell has 
a probability Pi to receive a ball. Then Aj (the expectation value of the number of the 
balls in the i-th. cell) reads Aj = Mpi with Y^iLi — 

The probability that the i-th cell receives rii balls is given by the binomial distribution. 
Its large M limit with Aj fixed is the Poisson distribution Px-{iT-i)'- 

BMA^i) = ,, Pr il-p^r-'''^P^An^) (M^oo). (A.2) 

Therefore, the probability that a certain image n = (ni, n2, ■ ■ ■, n^) is realized reads 

Pxin) = l[P^An^) = I[ (A.3) 

i=i i=i 

where the normalization is given by J2'^^=o P\iiP'i) = 1 (^ = 1; 2, ■ ■ ■ , A^). 

Since nj may be large as M is large, we introduce a small "quantum" q and define the 
finite image /j and the default model rrii as 

fi = qrii, rrii = q\i. (A.4) 

In terms of (|A.4|) , P(/ G V^) is written as 



r , df- ^ A"' e-^» r ^ df- e^^-^)/'' 

p(/ G V) = E ^a(^) - / n^^- / n%7f^. (A.5) 

•'^ ^ r=i nil JvfJ^Vfi (27rg)^/^ 

where small q is assumed for converting the sum to the integral, and the Stirling's formula, 
n! ~ \/27rne" ^°sn-n^ jg ^gg^j obtain the last expression. 
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S{f) is nothing but the Shannon- Jaynes entropy, 



s{f) = E 



N / r 

Ji 



fi-nii- fi log 

i=l . 

Comparing ( [A.5|) with ( |A.1| ), one finds 



(A.6) 



? = «-\ W] = I[% Zs{a)=r-^) . (A.7) 



Note also that the measure M{f) defined by 

N 

M{f) = n /r'/' (A.8) 

i=l 

is rewritten as 



M(/) = yd^ with = = (A.9) 

The metric tensor —Qij for the functional integral over / is thus related to the curvature 
matrix of S{f). 

B Axiomatic construction of entropy 

The axiomatic construction of the Shannon- Jaynes entropy S given below is a modified 



version of that given in |3|]. We have changed the statement of each axiom so that it 
becomes easier to understand the idea behind. 

For positive semi-definite distribution /(x), we want to assign a real number S{f) (the 
Shannon- Jaynes entropy) as 

/ is a more plausible image than g ^ 'S'(/) > S{g). (B.IO) 

If there exists an external constraint on f{x) such as C{f{x)) = 0, the most plausible 
image is given by the following condition 

6f[S{f)-XC{f)] = 0, (B.ll) 

with A being a Lagrange multiplier. The explicit form of S is uniquely fixed by the 
following axioms. 

Axiom I: Locality 

S{f) is a local functional of f{x) without derivatives. Namely, there is no correlation 
between the images at different x. 
This leads to a form 

S{f) = f dx m{x) e{f{x),x). (B.12) 



Here m(x) is a positive definite function which defines the integration measure. 6 is an 
arbitrary local function of f{x) and x without derivatives acting on /. 
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Axiom II: Coordinate Invariance 

f{x) and m(x) transform as scalar densities under tlie coordinate transformation x' 
x'{x), namely, f{x)dx = f\x')dx' and m{x)dx = m'{x')dx' . Also, S' is a scalar. 

This axiom allows only two invariants for constructing S in ( p.l2|) : m{x)dx 
m'{x')dx' and f{x)/m{x) = f'{x')/m'{x'). Therefore, 



S{f) = dx m{x) 0(/(x)/m(x)). 



(B.13) 



Axiom III: System Independence 

If X and y are two independent variables, the image F{x, y) is written as a product form 
F{x,y) = f{x)g{y) together with the integration measure m{x,y) = mf{x)mg{y). Fur- 
thermore, the first variation of S{F) with respect to F{x, y) leads to an additive form 
with some functions a{x) and P{y); 

6S{F) 



aix)+l3iy). 



(B.14) 



SF{x,y) 

From this axiom, images f{x) and g{y) are determined independently from the varia- 
tional equation 6S{F)/6F{x,y) = 0. First of all, ( [B.13|) reads 



S{F) = J dx J dy m{x,y) (j){F{x,y)/m{x,y)). (B.15) 

Acting the derivative d^/dxdy on eg . (P .141 ) with ( p.l5| ) and using Z = F{x, y)/m{x, y) = 
{f{x)/mf{x)){g{y)/mg{y)), one obtains 

where cr{Z) = d(j){Z)/dZ . The solution of this differential equation is <y{Z) = Ci logZ — Cq, 
which leads, up to an irrelevant constant, to 

(l){Z) =ciZ log Z -{co + ci)Z. (B.17) 

Thus one arrives at 

fix)' 



S{f) = dx m{x)(j){f /m) 



dx f{x) 



Ci log 



mix) 



(co + Ci) 



(B.18) 



where we have dropped the suffix / in mf{x) for simplicity. The curvature of 5* is dictated 
by ci, since {6/6fyS{f) = ci/ f and / > 0. We choose ci = —1 for having S bounded 
from above and for overall normalization. 



Axiom IV: Scaling 

If there is no external constraint on f{x), the initial measure is recovered after the varia- 
tion, i.e., /(x) = m(x). 

This axiom leads to cq = in ( p.l8| ), since the unconstrained solution of 6S{f)/6f = 
is f{x) = m(x)e'^°^^^. Also, it is convenient to add a constant to S{f) in ( p.l8|) to make 
S{f = m) = 0. Thus we arrive at 

fix) ^ 



Sif) 



dx 



fix) 



m(xi 



fix) log 



m(xi 



(B.19) 



which is the desired expression of the Shannon- Jaynes entropy. 
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C The singular value decomposition 

Here we give a proof of the singular value decomposition (SVD) of a general mxn matrix 



P following 1^ 



The singular values of a matrix P are defined as the square root of the eigenvalues 
of P^P. By definition, P^P is an n x n Hermitian matrix and has real and non-negative 
eigenvalues. We also define the norm of the vector x G C" and the spectral norm of the 
mxn matrix P, respectively, as 



X II2 
Ph 



nl/2 



U=l 

l^maximum eigenvalue of P^P 



1/2 



max{x'' P^ Px) ''' (x G C", || x ||2= 1). 



1/2 



(C.20) 

(C.21) 
(C.22) 



SVD Theorem: 

Let P be an m X matrix (m > n), f/ be an m x m unitary matrix, and ^ be an n x n 
unitary matrix. Then, P can be decomposed as 



P = UEV\ 



(C.23) 



where S is an m x n diagonal matrix with the diagonal elements being the singular values 

of P, namely, S = diag(^i, ^2, ■ ■ ■ , ^n) with ^1 > 6 > ■ • ■ > > 0. 

Proof: 

Define the maximum singular value of P as ^1. Then, there exists a vector xi which 
satisfies the relation: = (s|P''"Pxi). Also there exists a vector yi G C™(|| y \\2= 1) 
with the property, Pxi = ^lyi. 

For Xi and yi obtained above, one may introduce non-square matrices U2 and V2 such 
that Ui and Vi given below become anmxm unitary matrix and an n x tt, unitary matrix, 
respectively: 



U, = iy,,U2), V, = ix,,V2). 
Using Ui and Vi defined above, P is transformed into Pi as 



^1 4 
Q2 



(C.24) 



(C.25) 



where Zi G C" ^ and Q2 being an (m — 1) x (n — 1) matrix. 
Now, we show that zi is actually a null vector: 



P IP- 
II2" 



> 



^1 II 2 

T 



max X 



x) 



^1+ II ^1 112 
1 

ei+ II z. 111 

> e?+IUl|l2- 



II z. 



Q2Z1 1I2 



(C.26) 
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Since ,^1 is the maximum singular value of P and Pi, ^2 defined as the maximum SV of 
Q2 satisfies ^1 > ^2- Applying the same procedure to Q2, Q3, ■ ■ ■, one finds 



V 



J 



(C.27) 



Thus one arrives at the singular value decomposition P = UEV^ . (QED) 

One may neglect the irrelevant components of U and S so that they are an m x n 
matrix and an n x n matrix, respectively. In this case, U satisfies the condition U'^U = 1, 
while V^V = VV^ = 1. This form of SVD is used in the text. 



45 



References 



[1] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B147 (1979) 385, 
448. 

[2] M. A. Shifman, Prog. Theor. Phys. SuppL 131 (1998) 1. 

[3] T. Hatsuda and T. Kunihiro, Phys. Rev. Lett. 55 (1985) 158; Phys. Lett. B185 
(1987) 304. 

[4] M. Dey, V. L. Eletsky and B. L. loffe, Phys. Lett. B252 (1990) 620. 
C. Gale and J. Kapusta, NucL Phys. B357 (1991) 65. 
S. H. Lee, C. Song and H. Yabu, Phys. Lett. B341 (1995) 407. 

C. Song, P.W. Xia and C. M. Ko, Phys. Rev. C54 (1996) 3218. 

[5] M. Asakawa and C. M. Ko, Phys. Rev. C 48 (1993) R526. 

M. Herrmann, B. Friman and W. Norenberg, NucL Phys. A560 (1993) 411. 
B. Friman and H. J. Pirner, Nucl. Phys. A617 (1997) 496. 

W. Peters, M. Post, H. Lenske, S. Leupold and U. Mosel, Nucl. Phys. A632 (1998) 
109. 

G. E. Brown, G. Q. Li, R. Rapp, M. Rho and J. Wambach, Acta Phys. Pol. B29 
(1998) 2309. 

D. Cabrera, E. Oset and M. J. Vicente Vacas, |nucl-th/0011037 . 

[6] T. Hatsuda and T. Kunihiro, Phys. Rep. 247 (1994) 221. 
G. E. Brown and M. Rho, Phys. Rep. 269 (1996) 333. 

[7] E. V. Shuryak, Rev. Mod. Phys. 65 (1993) 1. 

[8] G. Agakichiev et al. (CERES Collaboration), Phys. Rev. Lett. 75 (1995) 1272; Phys. 
Lett. B422 (1998) 405. 

[9] M. C. Abreu et al. (NA50 Collaboration), Phys. Lett. B477 (2000) 28. 

[10] R. Rapp and J. Wambach, Adv. Nucl. Phys. 25 (2000) 1 (|hep-ph/9909"229| ). 

J. Alam, S. Sarkar, P. Roy, T. Hatsuda and B. Sinha, Ann. Phys. (N.Y.), 286 (2000) 

159 (|hep-ph/9909"26^ ). 

G. Chanfray, |nucl-th/0012"068| . 

[11] H. Satz, Rept. Prog. Phys. 63 (2000) 1511 ( |hep-ph/0007"069D . 

[12] Lattice 99 Proceedings, Nucl. Phys. B (Proc. Suppl.) 83 (2000) 1. 

[13] S. Aoki et al. (CP-PACS Collaboration), Phys. Rev. Lett. 84 (2000) 238. 

S. Aoki, Lattice Calculations and Hadron Physics, eConf C990809 (2000) 657 ( [hep- 
ph/9912"28g ). 



K. Kanaya, Hadronic Properties from Lattice QCD with Dynamical Quarks, [hep- 



ph/0005294 . 



46 



[14] M. -C. Chu, J. M. Grandy, S. Huang and J. W. Negele, Phys. Rev. D 48 (1993) 3340. 
D. B. Leinweber, Phys. Rev. D 51 (1995) 6369. 

D. Makovoz and G. A. Miller, Nucl. Phys. B468 (1996) 293. 
C. Allton and S. Capitani, Nucl. Phys. B526 (1998) 463. 

[15] T. Hashimoto, A. Nakamura and I. O. Stamatescu, Nucl. Phys. B400 (1993) 267; 
ibid. B406 (1993) 325. 

Ph. de Forcrand et al. (QCD-TARO Collaboration), |hep-lat/00080U5| . 

[16] Y. Nakahara, M. Asakawa and T. Hatsuda, Phys. Rev. D60 (Rapid Comm.) (1999) 
091503 ( |hep-lat/9905"034|) . 

Y. Nakahara M. Asakawa and T. Hatsuda, Nucl. Phys. (Proc. Suppl.) 83 (2000) 191 
( |hep-lat/9909137| ). 

[17] C. E. Shannon and W. Weaver, The Mathematical Theory of Communication, (Univ. 
of Illinois Press, Urbana, 1949). 

[18] E. T. Jaynes, Phys. Rev. 106 (1957) 620; ibid. 108 (1957) 171. 

[19] E.T. Jaynes, How does brain do plausible reasoning?, Stanford Univ. Microwave 
Lab. report 421 (1957), reprinted in Maximum- Entropy and Bayesian Methods in 
Science and Engineering, vol.1, pp. 1-24, eds. G. J. Erickson and C. R. Smith, (Kluwer 
Academic Publishers, London, 1988). 

[20] B. R. Frieden, J. Opt. Soc. Am., 62 (1972) 511. 

[21] N. Wu, The Maximum Entropy Method, (Springer- Verlag, Berlin, 1997). 

[22] See a recent review, M. Jarrell and J. E. Gubernatis, Phys. Rep. 269 (1996) 133. 

[23] R. N. Silver et al., Phys. Rev. Lett. 65 (1990) 496. 
R. N. Silver et al, Phys. Rev. B41 (1990) 2380. 
J. E. Gubernatis et al., Phys. Rev. B44 (1991) 6011. 

R. Preuss, W. Hanke and W. von der Linden, Phys. Rev. Lett. 75 (1995) 1344. 

W. von der Linden, R. Preuss and W. Hanke, J. Phys. 8 (1996) 3881. 

R. Preuss, W. Hanke, C. Grober and H. G. Evertz, Phys. Rev. Lett. 79 (1997) 1122. 

[24] S. E. Koonin, D. J. Dean and K. Langanke, Phys. Rep. 278 (1997) 1. 

[25] H. Jeffreys, Theory of Probability (Third Edition), (Oxford Univ. Press, Oxford, 
1998). 

G. E. P. Box and G. C. Tiao, Bayesian Inference in Statistical Analysis, (John Wiley 
and Sons, New York, 1992). 

[26] Ph. de Forcrand et al., Nucl. Phys. B (Proc. Suppl.) 63A-C (1998) 460. 

E. G. Klepfish, C. E. Creffield and E. P. Pike, Nucl. Phys. B (Proc. Suppl.) 63A-C 
(1998) 655. 

[27] J. Negele and H. Orland, Quantum Many-Particle Systems, (Addison- Wesley, New 
York, 1988). 



47 



[28] T. Hatsuda and S. H. Lee, Phys. Rev. C 46 (1992) R34. 

T. Hatsuda, Y. Koike and S. H. Lee, Nucl. Phys. B394, (1993) 221. 

S. H. Lee , Phys. Rev. C 57 (1998) 927; Erratum-ibid. C 58 (1998) 3771. 

There are some previous attempts to formulate the QCD sum rules at finite T by 

A. I. Bochkarev and M. E. Shaposhinokov, Nucl. Phys. B268 (1986) 220, and by 

T. Hashimoto, K. Hirose, T. Kanki and O. Miyamura, in the Proceedings of Santa 

Fe Strong Int. Workshop, (Santa Fe, April 8-11, 1986). In these works, however, key 

ingredients of the in-medium sum rules such as the existence of the tensor condensates 

at finite T and the consistent separation of the hard scale (M) and the soft scales (T 

and k) were not recognized. 

[29] M. Asakawa and C. M. Ko, Nucl. Phys. A572 (1994) 732. 

T. Hatsuda, S. H. Lee and H. Shiomi, Phys. Rev. C52 (1995) 3364. 

T. D. Cohen, R. J. Furnstahl, D. K. Griegel and X. Jin, Prog. Part. Nucl. Phys. 35 

(1995) 221. 

S. H. Lee, Nucl. Phys. A638 (1998) 183c, and references therein. 

S. Leupold and U. Mosel, Prog. Part. Nucl. Phys. 42 (1999) 221, and references 

therein. 

E. Marco and W. Weise, Phys. Lett. B482 (2000) 87, and references therein. 

[30] C. R. Smith and G. Erickson, in Maximum- Entropy and Bayesian Methods, pp. 29-44, 
ed. J. Skilling (Kluwer Academic Publishers, London, 1989). 

[31] J. Skilling, in Maximum Entropy and Bayesian Methods, pp. 45-52, ed. J. Skilling 
(Kluwer, Academic Pubhshers, London, 1989); S. F. Gull, ibid. pp. 53-71. 

[32] S. F. Gull and G. J. Daniell, Nature 272 (1978) 686. 

[33] E. T. Jaynes, in Maximum Entropy and Bayesian Methods in Applied Statistics, 
pp. 26-58, ed. J. H. Justice, (Cambridge Univ. Press, Cambridge, 1986). 

[34] J. Skilling, in Maximum Entropy and Bayesian Methods in Science and Engineering, 
vol.1, pp. 173-187, eds. G. J. Erickson and C. R. Smith, (Kluwer Academic Publishers, 
London, 1988). 

[35] R. K. Bryan, Eur. Biophys. J., 18 (1990) 165. 

[36] F. Chatelin, Eigenvalues of Matrices, (Wiley & Sons, New York, 1993). 

[37] W. H. Press et al.. Numerical Recipes, 2nd edition, (Cambridge Univ. Press, Cam- 
bridge, 1994). 

[38] For the comprehensive description on lattice gauge theories, see, for example, 

M. Creutz, Quarks, Gluons and Lattices, (Cambridge Univ. Press, Cambridge, 1983). 
I. Montvay and G. Miinster, Quantum Fields on a Lattice, (Cambridge Univ. Press, 
Cambridge, 1994). 

R. Gupta, Introduction to Lattice QCD, [hep-lat / 9807028| . 

M. Di Pierro, From Monte Carlo Integration to Lattice Quantum Chromo Dynamics, 
hep-iat/OOIMm . 



48 



[39] M. Creutz, Phys. Rev. D 21 (1980) 2308. 

A. D. Kennedy and B. J. Pendleton, Phys. Lett. B156 (1985) 393. 

F. R. Brown and T. J. Woch, Phys. Rev. Lett. 58 (1987) 2394. 
M. Creutz, Phys. Rev. D 36 (1987) 515. 

[40] The MILC code ver. 5, 

|http:/ /chodhna. cop. uop.edu/~hetrick/miIc . 

M. Asakawa, T. Hatsuda and Y. Nakahara, paper in preparation. 

M. Asakawa, T. Hatsuda, Y. Nakahara and S. Shibata, paper in preparation. 
S. Shibata, Master Thesis, Nagoya University (2000), in Japanese. 
See also, S. Sasaki, [hep-ph/0004254 

L.J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rep. 127 (1985) 1. 

M. G6ckeler et al., Nucl. Phys. B544 (1999) 699. 

D.E. Groom et al. (Particle Data Group), Eur. Phys. J. C15 (2000) 1. 

Y. Iwasaki et al., Phys. Rev. D 53 (1996) 6443. 
T. Bhattacharya et al., Phys. Rev. D 53 (1996) 6486. 

M. A. Shifman, A. I. Vainshtein and V. L Zakharov, Nucl. Phys. B147 (1979) 519. 
T. Hatsuda, E. M. Henley, T. Meissner and G. Krein, Phys. Rev. C49 (1994) 452. 
S-L. Zhu, W. Y. P. Hwang and Z-S. Yang, Phys. Rev. D57 (1998) 1524. 
S. Narison, Nucl. Phys. A675 (2000) 54. 

M. Asakawa, T. Hatsuda and Y. Nakahara, in progress. 

Proceedings of the workshop. Possible Existence of the a-Meson and its Implications 
to Hadron Physics, (June 12-14, 2000, Kyoto, Japan), eds. M. Ishida et al., KEK- 
proceedings/2000-4. 

See also, M. Alford and R. L. Jaffe, Nucl. Phys. B578 (2000) 367. 
R. D. Pisarski, Phys. Lett. BllO (1982) 155. 

G. E. Brown and M. Rho, Phys. Rev. Lett. 66 (1991) 2720. 

S. Chiku and T. Hatsuda, Phys. Rev. D58 (1998) 076001. 
T. Hatsuda, T. Kunihiro and H. Shimizu, Phys. Rev. Lett. 82 (1999) 2840. 
See also a short review, T. Hatsuda and T. Kunihiro, Chiral Symmetry Restoration 
and the a-meson, ( [hep-ph/0010039]) and references therein. 

T. Hashimoto, K. Hirose, T. Kanki and O. Miyamura, Phys. Rev. Lett. 57 (1986) 
2123. 

T. Matsui and H. Satz, Phys. Lett. B178 (1986) 416. 
M. Asakawa, T. Hatsuda and Y. Nakahara, in progress. 



49 



[56] M. Oevers, C. Davies and J. Shigemitsu, |hep-lat/ 0009031 



[57] C. DeTar, Phys. Rev. D32 (1985) 276. 

C. DeTar and J. B. Kogut, Phys. Rev. Lett. 59 (1987) 399; Phys. Rev. D36 (1987) 
2828. 



[58] A. Peshier and M. Thoma, Phys. Rev. Lett. 84 (2000) 841. 
F. Karsch, M. G. Mustafa and M. H. Thoma, |hep-ph/0007093 



L Wetzorke and F. Karsch, |hep-lat/0008"008 . 



50 



